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Abstract 

Models  for  mass  transport  inside  a  porous  SOFC  anode  were  developed  based  on  Fick’s  model  (FM),  the  dusty-gas  model  (DGM) 
and  the  Stefan-Maxwell  model  (SMM)  to  predict  the  concentration  overpotential.  All  models  were  validated  with  experimental  data  for 
F^-FFO-Ar  and  CO-CO2  systems.  The  effect  of  pore  size  on  all  model  predictions  was  discussed.  It  was  concluded  that  the  dusty-gas 
model  is  the  most  appropriate  model  to  simulate  gas  transport  phenomena  inside  a  SOFC  anode.  However,  this  model  requires  numerical 
solution,  whereas  Fick’s  and  Stefan-Maxwell’ s  do  not.  It  was  found  that  the  SMM,  rather  than  the  FM,  is  a  good  approximation  of  the 
dusty-gas  model  for  H2-H2O  system,  except  in  the  case  of  high  current  density,  low  H2  concentration  and  low  porosity,  where  only  the 
DGM  is  recommended.  For  the  CO-CO2  system,  there  is  no  simple  rule  for  selecting  an  alternate  model  to  DGM.  Depending  on  the  CO 
concentration,  porosity  and  current  density,  the  FM  or  the  SMM  could  be  used.  The  only  restriction  is  for  small  porosities  where  only  the 
DGM  should  be  used.  This  paper  also  demonstrated  that  only  the  DGM  is  recommended  for  a  multicomponent  system  (H2-H2O-CO-CO2). 
©  2003  Elsevier  Science  B.V.  All  rights  reserved. 
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1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  promise  electricity  gener¬ 
ating  processes  for  stationary  applications  in  the  mid-term 
future.  This  is  because  the  electrical  efficiency  typically 
achieved  in  a  SOFC  is  much  greater  than  that  obtained 
from  a  conventional  heat  engine  and  any  other  type  of  fuel 
cells  [1].  Furthermore,  additional  efficiency  can  be  gained 
by  adding  a  bottoming  cycle  to  recover  heat  from  hot  gas 
exhausted  from  a  SOFC.  Another  key  advantage  is  that  a 
SOFC  can  operate  with  CO,  which  is  a  poison  in  other  types 
of  fuel  cells.  Nowadays,  three  basic  designs  are  studied; 
they  are  cathode-,  electrolyte-  and  anode-supported  designs. 
The  latter  design  has  been  considered  to  be  better  than  any 
other  designs  because  it  can  operate  at  lower  temperature 
and  the  highest  performance  can  be  achieved  [2^4]. 

It  has  been  reported  [3,4]  that  the  performance  of  an 
anode-supported  SOFC  at  high  operating  current  density 
or  high  fuel  utilization  is  governed  by  the  concentration 
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overpotential  at  the  anode  side.  This  can  be  explained  by 
the  fact  that,  at  high  demands,  the  requirement  for  the 
reactant  exceeds  its  capability  to  diffuse  through  the  rel¬ 
atively  thick  porous  ceramic  anode  to  the  reaction  site  at 
the  vicinity  of  the  anode-electrolyte  interface  [5].  For  this 
reason,  the  accurate  prediction  of  the  anodic  concentration 
overpotential  is  extremely  important.  In  particular,  for  a 
tubular  anode-supported  SOFC  where  high  fuel  utilization 
is  required,  a  small  error  in  the  concentration  overpoten¬ 
tial  calculation  may  cause  a  dramatic  change  to  its  design 
performance. 

The  anodic  concentration  overpotential  considers  the  dif¬ 
ference  in  gas  concentrations  between  the  anode-electrolyte 
interface  and  the  bulk.  Mass  transport  models  inside  the 
porous  SOFC  anode  must  be  applied  to  estimate  gas  concen¬ 
trations  at  the  anode-electrolyte  interface.  However,  such  a 
mass  transport  model  needs  to  be  sophisticated  enough  so 
that  it  can  take  into  account  parameters  such  as  complex 
functions  of  temperature,  pressure,  gas  concentrations,  and 
the  physical  properties  of  SOFC  materials  like  porosity, 
tortuosity  and  pore  size  of  the  electrode  materials. 

The  transport  of  gaseous  components  through  porous  me¬ 
dia  has  been  extensively  studied  over  the  years  as  described 
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Nomenclature 

B0  pore  characteristics 

Def  overall  effective  diffusion  coefficient  of 
component  i  (cm2/s) 

D^im  molecular  diffusion  coefficient  of  component  i 
in  all  other  components  present  in  the  mixture 
(cm2/s) 

of  knudsen  diffusion  coefficient  of  component  i 
(cm2/s) 

DGM  dusty-gas  model 

F  Faraday’s  constant  (96500  C/mol) 

FM  Fick’s  model 

J  operating  current  density  (A/cm2) 

Mi  molecular  weight  of  component  i 

Ni  mol  flux  of  gas  species  i  (mol/(cm2  s)) 

P  operating  pressure  (atm) 

R  gas  constant  (8.31441  J/(molK)) 

r,  rate  of  reaction  inside  porous  media 

(mol/(m3  s)) 

SMM  Stefan-Maxwell  model 

T  operating  temperature  (K) 

yi  mole  fraction  of  component  i 

Greek  letters 

a  parameter  defined  by  Eq.  (18) 

e  porosity  (fraction) 

f]Conc  concentration  overpotential  (mV) 

/x  viscosity  (N/ms) 

r  tortuosity 

Subscripts  and  superscripts 

A/E  interface  between  anode  and  electrolyte 

B  binary  component  system 

T  ternary  component  system 


in  many  references  [6-8].  In  general,  mass  transport  of 
components  inside  porous  media  can  be  described  using  ei¬ 
ther  the  extended  Fick’s  model  (FM)  or  the  dusty-gas  model 
(DGM)  as  explained  in  [8].  Both  FM  and  the  DGM  are  mass 
transport  equations  taking  into  account  Knudsen  diffusion, 
molecular  diffusion  and  the  effect  of  a  finite  pressure  gra¬ 
dient.  The  key  feature  of  the  DGM  differs  from  that  of  FM 


in  that  the  flux  ratio  in  DGM  depends  on  the  square-root  of 
gas  molecular  weight.  From  a  literature  review,  although  the 
DGM  is  superior  to  FM  in  its  capability  to  predict  the  fluxes 
inside  porous  media  [9],  FM  is  more  frequently  used  be¬ 
cause  it  allows  explicit  analytical  expressions  to  be  derived 
for  fluxes.  Furthermore,  Veldsink  et  al.  [10]  claimed  that  the 
deviation  between  FM  and  the  DGM  can  be  neglected  for 
the  description  of  transport  and  reactions  in  a  porous  cata¬ 
lyst.  In  SOFC  modelling,  numerous  studies  on  gas  transport 
phenomena  [2,4,11,12]  have  been  carried  out  using  FM. 
Some  other  references  [13,14]  eliminated  out  the  effect 
of  Knudsen  diffusion  leaving  on  by  the  Stefan-Maxwell 
model  in  their  mass  transport  equations.  Only  Yakabe  et  al. 
[3]  and  Lehnert  et  al.  [5]  used  more  complex  mass  trans¬ 
port  models,  which  were  similar  to  the  DGM.  However, 
their  mass  transport  models  were  complex  and  required 
numerical  methods  to  obtain  the  solution.  Consequently, 
the  objective  of  this  paper  is  to  evaluate  the  various  mass 
transport  models  by  validating  them  with  the  measured  con¬ 
centration  overpotential  data  obtained  from  Yakabe  et  al. 
[3],  The  binary  systems  of  H2-H2O  and  CO-CO2  as  well 
as  the  ternary  systems  of  H2-H20-Ar  are  considered. 


2.  Model  development 

2. 1.  Modelling  assumptions 

The  modelling  assumptions  are  consistent  with  those  used 

by  Yakabe  et  al.  [3],  and  they  are  as  follows: 

1.  H2  and  CO  are  the  anodic  reactants.  Only  binary  gas 
mixtures  of  CO-CO2  and  H2-H2O  and  a  ternary  gas 
mixture  of  H2-H20-Ar  are  considered. 

2.  Gas  concentrations  are  uniform  along  the  fuel  channel. 
This  assumption  is  valid  because  the  length  of  the  fuel 
channel  is  very  small  (20  mm).  Fuel  utilization  is  less 
than  5%  and,  therefore,  assumed  to  be  constant.  Accord¬ 
ingly,  only  a  ID  model  in  the  z-direction  was  examined 
as  shown  in  Fig.  1. 

3.  The  current  density  is  uniform  over  the  entire  electrodes. 
This  assumption  is  valid  as  long  as  the  gas  concentration 
is  assumed  to  be  uniform  along  the  fuel  channel,  because 
the  current  density  produced  is  related  to  the  amount  of 
reactant  gas  diffusing  into  the  anode  layer. 

4.  Steady-state  condition. 


Fig.  1.  A  SOFC  geometry  used  in  model  development.  The  model  is  focused  on  diffusion  inside  a  porous  anode  layer. 
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5.  Temperature  and  total  pressure  are  constant  throughout 
the  entire  electrode.  The  total  pressure  gradient  inside  the 
electrode  is  assumed  to  be  negligible.  This  assumption  is 
quite  reasonable  since  there  is  no  net  change  in  the  num¬ 
ber  of  moles  in  the  gas  phase  [8]  due  to  electrochemical 
reactions.  Consider  the  electrochemical  reactions  of  H2 
and  CO  at  the  anode: 

H2  +  O2-  H20  +  2e“  (a) 

CO  +  O2-  C02  +  2e“  (b) 

6.  From  reactions  (a)  and  (b),  it  is  found  that  there  is  no  net 
change  in  the  number  of  moles  in  the  gas  phase. 

7.  Finally,  Lehnert  et  al.  [5]  reported  that  the  electrochem¬ 
ical  conversion  occurs  in  the  vicinity  of  the  electrolyte- 
anode  interface.  The  depth  of  the  reaction  zone  is  about 
50  p,m  above  the  anode-electrolyte  interface.  This  is  very 
small  compared  to  the  2  mm  anode  thickness  in  the  SOFC 
from  Yakabe  et  al.  [3j.  Therefore,  it  can  be  assumed  that 
the  electrochemical  reaction  is  a  heterogeneous  reaction, 
which  takes  place  at  the  anode-electrolyte  interface. 

2.2.  Mass  transport  equations 

Mass  transport  through  porous  medium  can  be  determined 
using  concepts  described  in  [8,10].  For  a  single  component, 
the  mass  transport  equation  can  be  written  as 

e  3(>'i  P) 

- =  —  V  ■  Ni  +  rj  (molm  3s  *)  (1) 

RT  dt 

where  £  and  Nj  represent  the  porosity  and  the  rate  of  mass 
transport,  respectively,  into  porous  media,  r,  is  the  rate  of  re¬ 
action  inside  the  porous  medium.  The  term  on  the  left-hand 
side  is  valid  when  unsteady  state  is  approached.  The  first 
and  second  terms  on  the  right-hand  side  represent  the  diffu¬ 
sion  rate  and  the  rate  of  reaction  inside  the  porous  medium. 
It  was  assumed  earlier  that  the  diffusion  process  is  at  steady 
state  and  that  the  electrochemical  reactions  take  place  at 
the  boundary  of  the  anode-electrolyte  interface  rather  than 
throughout  the  porous  medium.  Therefore,  within  the  an¬ 
ode,  only  the  first  term  on  the  right-hand  side  is  significant, 
Eq.  (1),  therefore,  becomes 

V  •  Ni  =  0  (2) 

The  rate  of  mass  transport,  Nj,  generally  depends  on  the 
operating  conditions  (reactant  concentration,  temperature 
and  pressure)  and  the  microstructure  of  material  (poro¬ 
sity,  tortuosity  and  pore  size).  Three  models  were  used  to 
develop  expressions  for  N,:  Fick’s,  the  dusty-gas  and  the 
Stefan-Maxwell  models. 

2.2.1.  Fick’s  model 

FM  is  the  simplest  form  used  to  describe  the  transport 
of  components  through  the  gas  phase  and  within  porous 
media.  The  general  extended  form  of  this  model  [10]  takes 
into  account  diffusion  and  convection  transport  and  is 


given  by: 


diyjP) 

dz 


BoVjP  dP\ 
F  3  z) 


(i  =  1, . . .  ,  n) 

(3) 


where  B0  is  the  permeability  coefficient,  Df(  the  effective 
diffusivity  of  species  i.  and  p.  the  viscosity  of  the  gas.  The 
first  and  second  terms  on  the  right-hand  side  account  for  dif¬ 
fusion  and  convection  transport,  respectively.  It  is  noted  that 
the  convection  transport  term  is  described  by  the  Darcy  equa¬ 
tion.  However,  it  is  assumed  that  the  total  pressure  change 
within  the  pore  is  insignificant  ( dP/dz  ~  0).  Thus,  only  the 
diffusion  transport  is  necessary  to  determine  the  rate  of  mass 
transport  Eq.  (3)  reduces  to: 


-of  3 (yiP) 

RT  dz 


(i  =  1,  . . .  ,  n) 


(4) 


The  diffusion  process  within  a  pore  can  be  typically  di¬ 
vided  into  two  diffusion  mechanisms:  molecular  diffusion 
and  Knudsen  diffusion.  Molecular  diffusion  is  dominant  for 
large  pore  sizes  and  high  system  pressures  while  Knudsen 
diffusion  becomes  significant  when  the  mean-free  path  of 
the  molecular  species  is  much  larger  than  the  pore  size. 
Therefore,  Df  can  be  written  by  combining  the  effective 
molecular  diffusion  (Df)  and  effective  Knudsen  diffusion 
(Of)  coefficients  as  follows  (Bosanquet  formula)  [10]: 


(5) 


Df  and  Df  depend  on  the  microstructure  of  the  porous 
anode  (porosity,  particle  size  and  tortuosity)  and  on  the  op¬ 
erating  conditions  (temperature  and  pressure).  The  calcula¬ 
tion  of  these  parameters  are  described  in  [3,4].  It  is  noted 
that  Df  given  in  Eq.  (5)  is  derived  from  the  assumption  of 
equimolar  counter  diffusion  of  reactant  and  product  species. 

For  binary  component  systems  (H2(l)-H20(2)  and 
C0(l)-C02(2)),  Eq.  (4)  becomes 


gtVdyi 

RT  dz 


(6) 


where  D'jllB  represents  the  effective  diffusion  for  binary  sys¬ 
tems  as  follows: 


Deff  - 

u\3  ~ 


1 


D' 


,eff 
1  ,k 


-1 


(7) 


At  the  anode-electrolyte  interface  (z  —  /a,  see  Fig.  1),  the 
amount  of  current  density  produced  is  governed  by  the  rate  of 
reactant  diffusing  into  the  porous  anode,  i.e.  N\  —  —J/2F, 
where  J  is  the  operating  current  density  and  F  the  Faraday’s 
constant.  Therefore, 


dyi  _  JRT 
d7  z=/a  “  ~2FDfPiP 


(8) 
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Inserting  Eq.  (6)  into  Eq.  (2)  and  solving  yields  to  the  fol¬ 
lowing  mole  fraction  profile: 

yi  =  Ci  +  Ciz  (9) 


where  C\  and  C2  are  constants.  By  using  the  Neuman  bound¬ 
ary  condition  at  z  =  /a  (Eq.  (8))  and  the  Dirichlet  boundary 
condition  y,-|z=o  =  yi.buik,  C 1  and  C2  are  determined  and 
Eq.  (9)  becomes: 


VI  =  yi.bulk  — 


JRT 

2  FDfBPZ 


(10) 


For  diffusion  with  heterogeneous  chemical  reaction,  the  flux 
ratios  are  governed  by  reaction  stoichiometry.  Summing 
Eq.  (14)  over  the  n  species  leads  to  the  Graham’s  law  of 
diffusion  in  gaseous  mixtures  [8]. 

n 

YJNifMi  =  0  (15) 

i— 1 

where  M-,  is  the  molecular  weight  of  component  i. 

For  two  binary  component  systems  (H2(l)-H20(2)  and 
C0(l)-C02(2)),  Eq.  (14)  becomes 


y 2  can  be  determined  by  the  fact  that  yi  +  v2  =  1 . 

For  ternary  component  systems  (H2(l)-H20(2)-Ar),  the 
formula  used  to  calculate  y\  and  y2  are  similar  to  the  one 
indicated  in  Eq.  (10)  except  that  the  overall  effective  diffu¬ 
sion  coefficient  for  the  ternary  component  system  (Z)|fj)  is 
derived  from  the  Stefan-Maxwell  relation.  The  mole  frac¬ 
tion  of  component  1  inside  the  porous  anode  is  expressed 
as  follows: 


vi  =  yi.bulk  — 


JRT 

2hT)fj  P  " 


(11) 


where 


neff  = _ i _ 

1,T  (l/£>f,)+(l/^iL)+((l/JDf)-(l/JDfl))(l-yAr) 

(12) 

and,  where  Dj®  is  the  effective  binary  diffusion  coefficient 
of  component  1  in  argon,  yAr  represents  the  mole  fraction 
of  argon  present  in  the  gaseous  mixtures. 


2.2.2.  Dusty -gas  model  (DGM) 

The  DGM  includes  the  Stefan-Maxwell  formulation  and 
takes  into  account  Knudsen  diffusion  [15].  It  is  assumed 
from  this  model  that  pore  walls  consist  of  giant  molecules 
(‘dust’)  uniformly  distributed  in  space.  These  dust  molecules 
are  considered  to  be  a  dummy,  or  pseudo,  species  in  the 
mixture.  The  general  form  of  the  DGM  is  shown  in  Eq.  (13): 


Ni 

Dfk 


E 


yjNi  -  yi  Nj 


Df 


1 

RT 


+  yi 


d  P 
d  z 


BoP  \ 


(13) 


The  second  term  on  the  right-hand  side  is  called  the  perme¬ 
ation  flux  and  is  applied  to  take  into  account  the  effect  of 
total  pressure  gradient  on  mass  transport.  However,  it  is  as¬ 
sumed  in  this  study  that  the  total  pressure  is  uniform  over 
the  entire  depth  of  the  porous  anode.  Thus,  the  definite  total 
pressure  gradient  term  can  be  ignored.  Consequently,  only 
the  diffusion  flux  is  addressed.  Eq.  (13)  reduces  to: 

Ni  1  y-'  yjNi  -  yiNi  p  n4. 

Def,f  \  .  Df  RT  dz  ’ 

;=1,  j¥=‘  V 


(Vi  yiN\  —  y\N2 
TF  T 


!)' 


,eff 


I) 


,eff 


1  ,k  ^12 

Rearranging  Eq.  (16)  gives, 


jP  dyi 
RT  dz 


N 1 


1 _ yz  -  yiNi/Ni^  _  P  dyi 


O' 


,eff 

IT 


D' 


,eff 

12 


RT  dz 


(16) 


(17) 


Because  y2  =  \  —  y\  and  N2/Ni  =  —^/M\/M2  (Graham’s 
law).  All  can  be  written  as  follow: 

-1 


P 

N 1  = - 

RT 


1  —  avi  1 


D' 


,eff 

12 


D 


,eff 

IT. 


dyi 

dz 


(18) 


where 

a  =  1  • 


Mi\ 

M2) 


1/2 


Substituting  Eq.  (18)  into  Eq.  (2)  gives 

-1-1 


dzyi 

dz2 


I) 


a 

]eff 

12 


1  —  ay  1  1 


iy 


eff 

12 


I) 


,eff 
1  ,k_ 


dyiV 

dz  / 


=  0 


(19) 


Eq.  (19)  is  in  the  form  of  an  ordinary  differential  equation, 
which  can  be  solved  by  using  the  two  following  initial  con¬ 
ditions  : 

IC.l  :  yiU=o  =  M.buik 


dvi 

JRT 

1 

—  ayi.buik 

1 

dz 

z=0 

2PF 

Deff 

^12 

neff 

ul,k  J 

Similarly,  for  ternary  component  systems  (H2(l)-H20(2)- 
Ar),  the  differential  equation  for  component  i  derived  from 
Eq.  (16)  is  as  follows 


d2yi 

dz2 


D' 


a 

(eff 

12 


D' 


1 

ieff 

IT 


D' 


ay  1 


I) 


1 

(eff 

1  Ar 
-1 


“  +  (1  —  jAr) 


dz  ) 


1 


1 


/yell  Tjeff 
,-^12  Ar  / 


(20) 


The  two  initial  conditions  are: 
IC.l  :  yi | -=o  =  A/ .bulk 
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IC.2  : 


d>7 

dz 


z=0 


JRT 

2PF 


1 


D 


eff 
1  ,k 


1 

— Ip  +  Cl-.VAr) 

u\  Ar 


1  \  «V1 

nclf  /  neff 
U\AiJ  U12 


2.3.1.  Fick’s  model 

The  concentration  overpotential  for  FM  can  be  addressed 
by  inserting  Eq.  (10)  for  the  binary  component  system  and 
Eq.  (11)  for  the  ternary  component  system  at  z  =  Za  into 
Eq.  (24): 


It  is  difficult  to  bnd  analytical  expressions  for  Eqs.  (19) 
and  (20).  Therefore,  the  differential  equations  were 
solved  numerically  in  MATLAB®  [16]  using  the  provided 
Runge-Kutta  method  ( ode45 ). 


„FM 
^Conc.B 


RT 

= - In 

2  F 


(yi,buik-(^T/2EDfBP)/a)(l-y1,buik)" 
yi,bulk(l+(TRr/2EDfBP)/a-yt,buik)  _ 

(25) 


2.2.3.  Stefan-Maxwell  model  (SMM) 

The  SMM  is  a  well-known  mass  transport  model  applied 
to  nonporous  medium.  When  the  Knudsen  diffusion  term  is 
neglected,  only  the  Stefan-Maxwell  diffusion  term  is  effec¬ 
tive.  Then,  Eq.  (14)  simplifies  to: 

”  yjNj  -  yiNj  P  dy,- 

^  £>eff  RT  dx  J 

j=i.jfi  y 

Eq.  (21)  can  be  solved  easily  using  the  same  concept  as  de¬ 
scribed  previously  for  FM.  In  this  case,  analytical  solutions 
for  binary  and  ternary  component  systems  can  be  obtained: 


y  i  =  yi.buik  — 


JRT 

2FIJfl  Zr 


(binary  component  system) 


(22) 


y  1  =  Jl.bulk  — 


/Pr(T>^r(l-yAr)  +  T>^,VAr). 

2  FPDfiDfto 


(ternary  component  system) 


(23) 


2.3.  Calculation  of  the  concentration  overpotential 


„FM 

'/Conc.T 


X  In 


RT 
~  ~2F 

(.Vi,  bulk  —  (/P7’/2PDjttTP)Za)(l  —  Vl,  bulk  —  VAr,  bulk) 

yi,buik(i+(^W<T^)Za-yi,buik-yAr,buik) 

(26) 


where  and  ZJ* ljl!|-  are  the  overall  effective  diffusion  co¬ 
efficients  for  binary  and  ternary  components,  respectively, 
as  defined  by  Eqs.  (7)  and  (12). 


2.3.2.  Dusty-gas  model 

For  the  dusty-gas  model,  Eq.  (19)  for  the  binary  compo¬ 
nent  system  and  Eq.  (20)  for  the  ternary  component  system 
were  solved  numerically  using  the  Runge-Kutta  method  to 
obtain  the  gas  compositions  at  the  reaction  site.  Afterwards, 
the  calculated  gas  compositions  were  inserted  into  Eq.  (24) 
to  calculate  the  concentration  overpotential. 


2.3.3.  Stefan-Maxwell  model 

The  concentration  overpotential  when  only  the  Stefan- 
Maxwell  diffusion  term  is  considered  is  determined  by  sub¬ 
stituting  Eqs.  (22)  and  (23)  into  Eq.  (24).  The  results  are  as 
follows: 


The  concentration  overpotential  used  to  validate  the  per¬ 
formance  of  various  mass  transport  models  is  the  difference 
between  ideal  and  real  cell  voltage  (when  corrected  for 
Ohmic  and  activation  overpotentials).  Using  the  Nernst 


SMM  „ 
yConc.B  — 


RT  f (yiMik-VRT/lFD^ P)Za)(l-yi,buik)" 
2 F  n  [  yi,buik(l+(7PP/2P^fP)Za-y1,bulk)  _ 

(27) 


SMM 
7  Conc.T 


RT]n  r  Or  bulk  -  (JRT(DfAiq  -  yAr)  +  P^yAr)/2PPDgz^r)Za)(l  -  vi,bulk  -  VAr.buik) 

2^  n  L  yi,bulk(l  +  {JRT(DfAr(  1  -  yAr)  +  DflyAr)/2FPDe"DfAr)f  -  yi.bulk  -  yAr.buik) 


(28) 


potential  and  mole  fractions  of  each  component  at  the  inter¬ 
faces  between  anode  and  electrolyte,  and  mole  fraction  of 
each  component  in  the  bulk,  the  concentration  overpotential 
is  defined  as: 


RT 

OCodc  =  —  tt;  In 

2  F 


yi,A/EV2.bulk 

.yi,bulky2,A/£. 


(24) 


where  yi.A/E  and  yi.A/E  are  mole  fractions  of  components 

1  and  2  at  the  interface  between  anode  and  electrolyte, 
respectively.  The  concentration  overpotential  for  binary  and 
ternary  components  for  each  model  is  calculated  by  insert¬ 
ing  the  mole  fractions  of  components  1  and  2  at  z  =  Za 
resulting  from  each  model  into  Eq.  (24). 


3.  Results  and  discussion 

3.1.  Model  validation 

All  input  parameters  used  in  this  model  calculation  were 
extracted  from  the  experimental  data  of  Yakabe  et  al.  [3]  as 
listed  in  Table  1. 

The  concentration  overpotentials  obtained  from  FM,  the 
DGM  and  the  SMM  are  compared  with  Yakabe’ s  exper¬ 
imental  data  as  illustrated  in  Figs.  2  and  3.  These  figures 
represent  experimental  data  and  calculated  concentration 
overpotential  for  the  Fb-FUO-Ar  and  CO-CCU  systems, 
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Table  1 


Input  parameters 


Parameter 

Symbol 

Value 

Operating  temperature  (°C) 

T 

750 

Operating  pressure(atm) 

P 

1 

Current  density  (A/cm2) 

J 

0.30,  0.70,  1.0  for 
FF-FFO~Ar  system 
0.10,  0.30,  0.50  for 
CO-CO2  system 

Porosity  (fraction) 

£ 

0.46 

Average  pore  size  (p,m) 

d 

2.6 

Tortuosity 

T 

4.5  (fitted  to 
experimental  data) 

Anode  thickness  (mm) 

4 

2 

respectively.  The  aim  of  this  experiment  is  to  investigate 
the  effect  of  reactant  concentration  on  the  concentration 
overpotential.  The  vertical  axis  represents  the  concentra¬ 
tion  overpotential  (ijconc)  minus  its  basis  concentration 


overpotential  (/? conco)-  The  horizontal  axis  is  the  reactant 
concentration.  According  to  Yakabe  et  al.  [3],  f/conc()  is 
measured  at  H2/(H2  +  H20  +  Ar)  =  0.8  for  H2-H20-Ar 
and  at  CO/ (CO  +  C02)  =  0.64  for  C0-C02  system.  For 
H2-H20-Ar  system,  a  ratio  of  H2:H20  is  kept  at  4: 1  and  the 
concentration  of  H2  in  the  system  is  varied  by  the  degree  of 
dilution  of  Argon.  This  means  that  ljconcO  is  measured  when 
no  argon  is  present  in  the  system.  In  this  simulation,  the 
tortuosity  factor  is  treated  as  a  variable  parameter  necessary 
to  fit  the  experimental  results.  The  results  from  H2-H?0-Ar 
system  as  illustrated  in  Fig.  2  are  used  to  identify  the  tor¬ 
tuosity  factor.  The  best  fitted  tortuosity  factor  is  about  4.5, 
which  is  in  the  range  of  reported  value  between  3  and  6  [13] 
and  is  the  same  as  the  value  reported  by  Yakabe  et  al.[3]. 
From  Figs.  2  and  3,  f/conc  increases  as  concentrations  of 
H2  and  CO  decrease.  The  increase  of  rjconc  for  low  reactant 
concentration  is  greater  than  that  for  high  reactant  concen¬ 
tration.  In  other  words,  the  change  in  »;conc  for  low  reactant 
concentration  is  more  sensitive  to  the  change  of  reactant 


Fig.  2.  Comparison  between  the  experimental  data  and  calculated  concentration  overpotential  for  the  FF-FFO-Ar  system:  symbols  represent  experimental 
data,  lines  simulated  ones;  (♦)  0.3 A/cm2,  (A)  0.7 A/cm2,  (■)  1.0 A/cm2;  ( — )  FM,  ( - )  DGM,  ( — )  SMM. 


0.1  0.2  0.3  0.4  0.5  0.6  0.7 

C0/(C0+C02) 


Fig.  3.  Comparison  of  the  experimented  and  calculated  concentration  overpotential  for  the  CO-CO2  system:  symbols  represent  experimental  data,  lines 
simulated  ones;  (♦)  0.1  A/cm2,  (A)  0.3  A/cm2,  (■)  0.5  A/cm2;  ( — )  FM,  (-■■-)  DGM,  ( — )  SMM. 
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concentration  than  that  for  high  reactant  concentration. 
When  compared  among  current  densities,  the  increase  in 
current  density  results  in  increasing  ?;conc-  The  explanation 
is  that  f]conc  has  a  greater  impact  at  low  fuel  concentra¬ 
tions  and  high  operating  current  densities  than  at  high  fuel 
concentrations  and  low  operating  current  densities.  This  is 
because  at  high  fuel  concentration  and  low  operating  current 
densities,  the  bulk  gas  contains  sufficient  reactants  to  diffuse 
directly  to  the  anode-electrolyte  interface.  Consequently, 
the  amount  of  reactant  gas  at  the  anode-electrolyte  interface 
is  not  much  lower  than  that  in  the  bulk  gas.  This  results  in 
lowering  pconc  •  On  the  other  hand,  at  low  fuel  concentration 
and  high  operating  current  density,  the  amount  of  reactant 
in  the  bulk  gas  is  not  sufficient  to  attain  the  desired  operat¬ 
ing  current  density.  Then,  the  reactants  diffuse  slowly  to  the 
anode-electrolyte  interface.  This  results  in  starving  reac¬ 
tants  at  the  reaction  site  (anode-electrolyte  interface).  It  can 
be  noted  from  experimental  and  model  results  that  there  is  a 
low  concentration  limit  depending  on  the  operating  current 
density  below  which  r/conc  increases  dramatically.  When 
concentrations  of  kb  and  CO  are  near  this  threshold  limit,  a 
small  decrease  in  gas  concentration  has  a  tremendous  effect 
on  the  ?7conc-  This  means  that  the  accuracy  in  rj conc  pre¬ 
diction  is  crucial  for  low  reactant  concentration.  Therefore, 
although  the  fuel  utilization  from  Yakabe’s  experiment  is 
less  than  5%,  the  assumption  of  constant  current  density  is 
not  valid  when  the  fuel  concentration  is  close  to  its  thresh¬ 
old  limit.  Therefore,  the  1-D  model  is  adequate  only  for 
fuel  concentrations  much  greater  than  the  threshold  limit. 
However,  near  the  threshold  limit  a  2-D  model  is  required. 

From  Fig.  2,  it  is  found  that  for  H2-H20-Ar  system,  the 
performance  of  the  DGM  in  predicting  experimented  r;conc 
is  as  good  as  that  of  SMM  for  high  and  moderate  H2  concen¬ 
tration.  However,  its  performance  is  superior  to  that  of  the 
SMM  for  low  H2  concentration.  FM,  on  the  other  hand,  pro¬ 
vides  the  worst  matched  results  especially  at  low  H2  concen¬ 
trations.  The  reason  is  that  FM  uses  the  concept  of  equimo¬ 
lar  counter  diffusion  (  flux  ratio  is  equal  to  —1)  to  derive  the 
overall  effective  diffusion  coefficient  as  indicated  in  Eq.  (5). 
This  concept  is  no  longer  valid  when  Knudsen  diffusion  is 
dominant.  In  this  case,  Graham’s  law,  which  is  derived  from 
the  DGM,  must  be  applied  to  calculate  flux  ratios.  Graham’s 
law  states  that  the  flux  ratio  of  diffusion  inside  porous  media 
relies  inversely  on  the  square-root  of  molecular  weight  of  gas 
at  which  the  flux  ratio  is  calculated.  For  the  H2-H20-Ar  sys¬ 
tem,  (Vht/-Yh20  =  —  V18/2  =  —3  (argon  is  a  stagnant  gas), 
which  is  very  different  from  —1.  Therefore,  the  assump¬ 
tion  of  equimolar  counter  diffusion  is  not  valid.  The  SMM, 
on  the  other  hand,  can  satisfactory  assume  that  equimolar 
counter  diffusion  occurs  since  the  porous  medium  is  treated 
as  a  stagnant  gas  film.  The  Knudsen  diffusion  term  is  not 
taken  into  account  in  the  SMM.  However,  this  model  is  not 
applicable  for  modeling  gas  transport  phenomena  through 
small  pore  sizes.  Fortunately,  the  effect  of  Knudsen  diffu¬ 
sion  from  this  model  validation  is  not  so  critical  when  com¬ 
pared  to  molecular  diffusion  because  the  pore  size  is  so  large 


that  the  Knudsen  diffusivity  (Dh, ,k  ^  28.5cm2/s)  is  about 
four  times  greater  than  the  binary  diffusivity  (£>h,-h2o  ~ 
7.1cm2/s).  However,  the  molecular  diffusion  depends  not 
only  on  binary  diffusivity  but  reactant  concentration  as  well. 
The  low  molecular  diffusion  effect  is  obtained  for  low  reac¬ 
tant  concentration.  Therefore,  the  Knudsen  diffusion  can  be 
crucial  to  the  overall  diffusion  for  low  reactant  concentra¬ 
tion.  This  explains  the  reason  why  the  DGM  is  better  than 
the  SMM  for  predicting  rjconc  at  low  H2  concentrations. 

For  the  CO-CO2  system,  as  illustrated  in  Fig.  3,  the  per¬ 
formance  of  FM  and  the  SMM  is  as  good  as  that  of  the 
DGM  over  the  whole  range  of  CO  concentration  at  low 
current  density  (0.1  A/cm2).  For  moderate  current  densities 
(0.3  A/cm2),  both  the  DGM  and  FM  are  somewhat  better 
than  that  of  the  SMM  especially  at  low  concentration.  This 
can  be  explained  by  the  fact  that  the  flux  ratio  of  CO  to 
CO2  is  closer  to  —1  (Nco/Nco2  =  —  V44/2 8  =  —1.25). 
This  means  that  the  assumption  of  equimolar  counter  diffu¬ 
sion  used  in  FM  becomes  more  valid.  In  addition,  FM  takes 
into  account  Knudsen  diffusion,  which  is  not  considered  for 
SMM.  In  Fig.  3,  it  is  also  shown  that  all  models  provide 
bad  predictions  of  the  measured  data  at  a  current  density 
of  0.50  A/cm2.  One  possible  reason  for  this  is  that  at  high 
current  density,  the  concentration  profile  of  CO  along  fuel 
channel  is  so  different  that  the  assumption  of  uniform  gas 
distribution  is  in  valid  since  the  diffusion  rate  of  CO  is  slow. 

3.2.  Effect  of  pore  size 

The  influence  of  pore  size  is  one  of  the  key  parameters 
used  to  evaluate  different  mass  transport  models.  This  is 
because  an  appropriate  mass  transport  model  for  large  pore 
size  may  be  not  suitable  for  small  pore  size.  Figs.  4  and 
5  depict  the  relationship  between  pore  size  and  the  con¬ 
centration  overpotential  (rjconc)  for  H2-H2O  and  CO-CO2 
systems,  respectively.  The  current  densities  remain  at 
1  A/cm2  for  H2-H2O  system  and  0.30  A/cm2  for  CO-CO2 
system.  The  typical  range  of  pore  size  is  varied  from  1  p,m 
(standard  size  of  anode  material,  [5])  to  2.6  pm  [3].  This 
investigation  is  performed  at  high,  moderate  and  low  reac¬ 
tant  concentrations.  The  mole  fractions  of  reactants  are  kept 
at  0.8,  0.4  and  0.25  for  high,  moderate  and  low  reactant 
concentrations,  respectively.  In  Figs.  4  and  5,  it  is  obvious 
that  pore  size  has  no  effect  with  using  the  SMM,  which  is 
to  be  expected  since  Knudsen  diffusion  is  not  taken  into 
consideration.  The  variable  rjconc  estimated  from  both  FM 
and  the  DGM,  in  contrast,  have  the  same  pattern  where  it 
increases  as  pore  size  decreases.  This  means  that  when  pore 
size  is  reduced,  the  Knudsen  diffusion  becomes  predomi¬ 
nant  due  to  the  decreasing  capability  of  gas  to  diffuse  inside 
the  porous  structure.  Therefore,  the  r/ccmc  estimated  from 
both  FM  and  DGM  should  be  equal  to  or  higher  than  the 
?7Conc  obtained  from  the  SMM.  However,  it  is  observed  that 
fjConc  calculated  from  the  DGM  is  below  >/conc  obtained 
from  the  SMM  for  large  pore  size  and  high  H2  concentra¬ 
tion  as  shown  in  Fig.  4(a)  in  the  case  of  H2-H2O  system. 
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Fig.  4.  Relationship  between  pore  size  and  the  concentration  overpo¬ 
tential  (t) cone)  for  the  H2-H2O  system  at  current  density  =  1  A/cm2, 
(a)  H2/(H2  +  H20)  =  0.80;  (b)  H2/(H2  +  H20)  =  0.40;  (c) 
H2/(H2  +  H20)  =  0.25. 
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Fig.  5.  Relationship  between  pore  size  and  the  concentration  overpoten¬ 
tial  (fjeone)  for  the  CO-CO2  system  at  current  density  =  0.30A/cm2. 
(a)  CO/(CO  +  C02)  =  0.80;  (b)  CO/(CO  +  C02)  =  0.40;  (c) 
CO/(CO  +  C02)  =  0.25. 


To  explain  Fig.  4(a)  consider  the  expressions  of  the  flux  for 
both  the  DGM  and  the  SMM  models: 


Nh2 


RT  \  Deff 

\  UH2,k 

for  the  DGM 


(1-0.67.vh2)\ 

Deff  / 

'ih-IFO  / 


d>’H2 

dz 


(29) 


It  is  observed  that  there  are  two  parameters  that  make  the 
DGM  different  from  the  SMM;  k  and  vi  p .  For  large  pore 
sizes  and  high  concentrations,  Knudsen  diffusion  becomes 
negligible  (1/D^  k  -a-  0)  and  Eq.  (29)  becomes: 


Nh2 


■ph2-h2o^>  dyH2 
RT(l  —  0.67vh2)  d z 


for  the  DGM 


(31) 


Nh2 


D 


eff 

h2-h2o 


pdyH2 


dz 


for  the  SMM 


Comparing  Eqs.  (30)  and  (31),  it  is  apparent  that  the  hydro¬ 
gen  molar  flux  is  greater  when  determined  with  the  DGM 


RT 
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Table  2 

Criteria  for  choosing  appropriate  mass  transport  models  for  CO-CO2  system 


CO  concentration  (%)  Low  current  density  (0.1  A/cm2)  High  current  density  (0.3  A/cm2) 

Pore  size  Pore  size 

Small  (<1.5  p,m)  Intermediate  (1.5-2. 0  p,m)  Large  (>2.0  p.m)  Small  (<1.5  p-m)  Intermediate  ( 1 .5-2.0  pm)  Large  (>2.0  pm) 


Low  (25%) 

FM 

FM 

FM 

DGM  only 

DGM  only 

FM 

Intermediate  (40%) 

FM 

FM 

FM 

FM 

FM 

FM 

High  (80%) 

DGM  only 

SMM 

SMM 

DGM  only 

SMM 

SMM 

Remark'.  DGM  is  the  most  accurate  model  in  every  case.  When  FM  or  SMM  is  indicated  in  the  table,  it  means  that  the  results  from  these  models  are 
within  15%  of  the  results  obtained  from  DGM. 


than  with  the  SMM.  Therefore,  for  large  pore  size,  r/conc  is 
lower  for  the  DGM  than  for  the  SMM. 

This  behaviour  does  not  occur  in  the  CO-CCD  system 
as  can  be  observed  in  Fig.  5.  This  is  because  the  flux  ratio 
used  in  the  DGM  is  almost  similar  to  that  used  in  SMM. 
Furthermore,  the  CO  molecule  is  much  larger  than  that  of 
FF.  Therefore,  Knudsen  diffusion  has  more  impact  on  the 
overall  diffusion  rate  even  at  large  pore  size.  It  is  also  noted 
that  FM  always  predicts  >)conc  higher  than  the  DGM.  The 
reason  for  this  is  that  the  flux  ratio  used  in  FM  is  lower 
than  that  used  in  DGM. 

From  Fig.  4(a)  and  (b),  it  is  shown  that  discrepancies 
between  DGM  and  SMM  are  smaller  than  that  between 
the  DGM  and  FM  over  the  entire  range  of  pore  size.  This 
means  that  the  SMM  is  a  good  approximation  of  the  DGM 
for  high  and  moderate  FF  concentrations.  At  low  H2  con¬ 
centrations,  as  exhibited  in  Fig.  4(c),  it  was  found  that  the 
SMM  is  an  appropriate  model  for  only  moderate  and  large 
pore  size.  This  is  not  an  appropriate  model  for  small  pore 
size  since  the  impact  of  Knudsen  diffusion  overcomes  the 
molecular  diffusion.  In  this  case,  it  is  recommended  to  use 
the  DGM  rather  than  either  FM  or  the  SMM.  At  high  CO 
concentrations  in  the  CO-CO2  system,  as  illustrated  in 
Fig.  5(a),  the  performance  of  the  SMM  in  predicting  i?conc 
is  superior  to  that  of  FM  for  moderate  and  large  pore  sizes. 
On  the  other  hand,  its  performance  is  a  somewhat  lower 
than  the  performance  of  FM  for  small  pore  sizes.  However, 
either  FM  or  the  SMM  can  be  applied  without  significant 
loss  of  accuracy.  At  moderate  CO  concentrations,  as  shown 
in  Fig.  5(b),  it  seems  that  the  difference  between  the  DGM 
and  FM  is  smaller  than  that  between  DGM  and  SMM  over 
the  whole  range  of  pore  size.  Consequently,  it  is  recom¬ 
mended  to  use  FM  in  moderate  CO  concentrations.  At  low 
CO  concentrations,  as  depicted  in  Fig.  5(c),  FM  is  a  good 
approximation  of  the  DGM  for  large  pore  size.  However, 
as  pore  size  decreases,  both  FM  and  the  SMM  show  large 
deviations  compared  to  the  DGM.  Therefore,  it  is  recom¬ 
mended  to  use  the  DGM  for  small  and  moderate  pore  sizes. 

3.3.  Criteria  for  choosing  an  appropriate 
mass  transport  model 

From  the  performance  comparison  of  the  DGM,  FM  and 
the  SMM  for  predicting  )]Conc,  it  is  found  that  there  are 


three  key  parameters  specifying  the  criteria  for  choosing  an 
appropriate  mass  transport  model:  current  density,  reactant 
concentration  and  pore  size.  Although  the  DGM  can  predict 
i?Conc  accurately,  it  requires  a  numerical  method  to  solve. 
Based  on  this  study,  it  can  be  concluded  that  SMM  is  a  good 
approximation  of  the  DGM  for  the  H2-H2O  system  except 
in  the  case  of  high  current  densities,  low  hydrogen  concen¬ 
trations  and  small  pore  sizes  where  only  the  DGM  is  recom¬ 
mended.  For  the  CO-CO2  system,  the  criteria  for  choosing 
the  appropriate  mass  transport  models  are  summarized  in 
Table  2. 

For  multicomponent  systems  such  as  the  H2-H2O-CO- 
CO2  system,  the  flux  ratios  among  gas  species  are  so  compli¬ 
cated  that  the  equimolar  counter  diffusion  assumption  is  no 
longer  valid.  In  this  case,  it  is  recommended  to  use  the  DGM 
rather  than  either  FM  or  the  SMM.  According  to  the  DGM, 
as  defined  in  Eq.  ( 14),  to  derive  the  flux  of  each  gas  species  in 
H2-H2O-CO-CO2  system,  it  is  needed  to  determine  flux  ra¬ 
tios  of  Nh2/Nyi2o,  Nco/Nco2’  Nco/Nn2’  NCo/Nh2o,  and 
NcoJNh2.  From  Graham’s  law  of  diffusion  /V|[2/;V|[,o  = 
—3  and  Nco/Nqo2  =  —1-25.  The  electrochemical  reactions 
of  H2  (H2  +  02“  -*  H20  +  2e“)  and  CO  (CO  +  O2”  -* 
CO2  +  2e_)  occur  simultaneously.  Therefore,  the  flux  ratio 
of  CO  to  Ht  can  be  related  to  the  current  density  produced 
by  CO  (Jco)  that  generated  by  H2  (./[[,  ),  i.e.  (Vco /  A')  [,  = 
Jco/Jn2-  Finally,  NCo/Nh2o  and  /Vco2/A'h2  are  functions 
of  both  Jco  and  -Al  and  are  determined  to  be  —  3Jco/-Ab 
and  —  0.8Jco/-Al-  respectively.  From  this  flux  ratio  calcula¬ 
tion,  it  is  demonstrated  that  the  equimolar  counter  diffusion 
approach  cannot  be  applied  to  this  system. 


4.  Conclusions 

Mass  transport  models  based  on  FM,  the  DGM  and  the 
SMM  were  developed  to  predict  the  concentration  over¬ 
potential  inside  a  SOFC  anode.  There  are  a  variety  of 
differences  between  the  models.  First  of  all,  the  SMM  does 
not  include  the  Knudsen  diffusion  term  to  account  for  pore 
size  effect  while  both  FM  and  the  DGM  do.  Secondly,  the 
equimolar  counter  diffusion  is  assumed  to  determine  the 
flux  ratio  in  both  FM  and  the  SMM.  By  contrast,  the  flux 
ratio  in  the  dusty-gas  model  is  calculated  from  the  ratio  of 
the  square-root  of  the  gas  molecular  weight.  This  is  derived 
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from  Graham’s  law  of  diffusion.  Finally,  FM  and  the  SMM 
can  be  derived  analytically  while  the  DGM  requires  a 
numerical  solution. 

It  was  found  that  three  parameters  define  the  appropriate 
mass  transport  model:  current  density,  reactant  concentration 
and  pore  size.  The  DGM  is  the  most  suitable  model  for  the 
H2-H2O  and  CO-CO2  systems.  This  is  because  it  takes  into 
account  Knudsen  diffusion  effect  as  well  as  Graham’s  law 
of  diffusion  to  calculate  the  flux  ratios.  However,  because 
of  its  complexity,  this  model  is  only  required  for  conditions 
of  high  operating  current  density,  low  reactant  concentration 
and  small  pore  size  where  high  accuracy  of  model  prediction 
is  required,  otherwise,  the  SMM  is  a  satisfactorily  appropri¬ 
ate  model  for  the  H2-H2O  system  while  the  FM  is  suitable 
for  the  CO-CO2  system.  For  the  multicomponent  system 
of  H2-H2O-CO-CO2;  however,  the  DGM  is  recommended. 
This  is  because  the  flux  ratio  calculation  is  complicated  and 
equimolar  counter  diffusion  cannot  be  assumed.  Further¬ 
more,  the  Knudsen  diffusion  effect  must  be  considered. 
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